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ABSTRACT 

We study the stability properties of Rayleigh unstable flows both in the purely hy¬ 
drodynamic and magnetohydrodynamic (MHD) regimes for two different values of 
the shear q - 2.1,4.2 (q — -dlnCl/dlnr) and compare it with the Keplerian case 
q — 1.5. We find that the q > 2 regime is unstable both in the hydrodynamic and in 
the MHD limit (with an initially weak magnetic field). In this regime, the velocity 
fluctuations dominate the magnetic fluctuations. In contrast, in the q <2 (magnetoro- 
tational instability (MRI)) regime the magnetic fluctuations dominate. This highlights 
two different paths to MHD turbulence implied by the two regimes, suggesting that 
in the q > 2 regime the instability produces primarily velocity fluctuations that cause 
magnetic fluctuations, with the causality reversed for the q <2 MRI unstable regime. 
We also find that the magnetic field correlation is increasingly localized as the shear is 
increased in the Rayleigh unstable regime. In calculating the time evolution of spatial 
averages of different terms in the MHD equations, we find that the q > 2 regime is 
dominated by terms which are nonlinear in the fluctuations, whereas for q < 2, the 
linear terms play a more significant role. 
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1 INTRODUCTION 

Differentially rotating flows are ubiquitous in astrophysics 
and studying their stability has been a long-standing 
enterprise. Using the local shearing box approximation 
(Goldreich & Lynden-Bell (1965), Hawley et al. (1995) with 
Keplerian shear (q = 1.5), numerical simulations have shown 
that the Magnetorotational Instability (MRI) leads to tur¬ 
bulent growth of stresses in the presence of a weak mag¬ 
netic field (for example, Velikhov (1959), Chandrasekhar 
(1960), Balbus & Hawley (1991)). The Rayleigh criterion, 
based on a linear modal analysis of axisymmetric pertur¬ 
bations, suggests that Keplerian flow is stable in hydro¬ 
dynamics. This, however, does not rule out the possibility 
of subcritical transition to turbulence (Balbus et al. (1996), 
Lesur & Longaretti (2005)). 

The (Rayleigh stable) Keplerian flow has understand¬ 
ably received the most attention because of its direct appli¬ 
cation in accretion discs, but here we focus on the stabil¬ 
ity properties of hydrodynamic and magnetohydrodynamic 
(MHD) flow in the Rayleigh unstable regime q > 2. A study 
of the Rayleigh unstable regime is of interest because a com¬ 
prehensive understanding of shear driven MHD turbulence 


requires knowing the differences in the q < 2 and q > 2 
regimes. Additionally, certain astrophysical flows are actu¬ 
ally thought to be Rayleigh unstable. These include counter 
rotating accretion discs (e.g., Dyda et al. (2015)), counter 
rotating galaxies (e.g., Corsini (2014)) and the plunging re¬ 
gion close to a black hole (e.g., Abramowicz et al. (1978), 
Abramowicz et al. (1996), Gammie (2004), Balbus (2012), 
Penna et al. (2013)). 

While the standard shearing box in the Rayleigh unsta¬ 
ble regime poses challenges that we discuss further in section 
2.2, certain properties of the shear instabilities in both the 
hydrodynamic and magnetohydrodynamic (MHD) case can 
be studied numerically with an appropriate configuration 
and code. Toward this end, we have conducted numerical 
simulations for three different values of q (1.5,2.1,4.2) both 
in pure hydrodynamics and MHD. We first used the publicly 
available finite volume code ATHENA 1 ((Gardiner & Stone 
(2005), Stone et al. (2008), Stone & Gardiner (2010)) and 
found that even though we started out with zero initial mo¬ 
menta, truncation errors introduced perturbations that led 
to the exponential growth of the mean momentum and the 
eventual crash of the simulation (the time step is inversely 
proportional to maximum velocity). We then chose the 
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pseudospectral code snoopy 2 (Lesur & Longaretti (2005), 
Lesur & Longaretti (2011)) to simulate q > 2, which con¬ 
serves the k = 0 mode. 

In section 2 we review the linear stability theory of 
hydrodynamic and magnetohydrodynamic shear flows and 
discuss it in the context of shearing box approximation. In 
section 3, we describe the numerical setup and simulation 
results. We conclude in section 4. 


2 STABILITY OF SHEAR FLOWS 


2.1 Linear analysis 

Following the discussion in Shakura & Postnov (2015), 
the dispersion relation for local axisymmetric perturba¬ 
tions of the form e' ( ' at ~ krr ~ kzi> (see also Balbus (2012) for 
the special case of k r = 0) with the initial magnetic 
field Bq pointing in the z direction is (Velikhov (1959), 
Chandrasekhar (1960), Balbus & Hawley (1991), Kato et al. 
(1998), Shakura & Postnov (2015)): 



(1) 


where k 2 = k 2 + k 2 , k 2 = 4£! 2 + rdQ. 2 /dr = 2Q 2 (2 - q), v 2 = 
Bq/(4ttpq) and po = initial density. The solution is 
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For the classical Rayleigh criterion in hydrodynamics, v A = 0 
and the above relation gives a> 2 = ( k z /k) 2 K 2 . This implies that 
purely hydrodynamic perturbations are stable as long as k 2 > 
0, or equivalently q <2. However, the addition of magnetic 
fields makes the q <2 regime unstable and instead a»^[ RI ~ 
(k 2 v A )/(K 2 dkl 2 /din r) = -q/(2 - q)k 2 v^ in the limit lrv\ « 1 
(Balbus 2012). 

We focus our attention to the q > 2 or k 2 < 0 regime 
in this paper. It is convenient to define the two different 
branches of Eq. 2 in the limit of k 2 v 2 A « 1 as: 
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where cjr = Rayleigh mode and u>vc = 
Velikhov-Chandrasekhar mode. As explained by 
Shakura & Postnov (2015), these modes are so named 
because we recover the classical Rayleigh instability crite¬ 
rion from the Rayleigh mode in the absence of magnetic 
field (va = 0), and the VC mode vanishes in this limit. In 
the regime tc < 0, it follows from above that the VC mode 
is stable for all wavenumbers and only the Rayleigh mode 
is unstable. This distinction between the Rayleigh and VC 
mode was not made in Balbus (2012). 


2 http://ipag.osug.fr/~lesurg/snoopy.html 


2.2 Shearing box in the Rayleigh unstable regime 


The shearing box approximation in the ideal compressible 
MHD limit is discussed in Nauman & Blackman (2015). 
Here we revisit that discussion in the context of non-ideal 
incompressible MHD equations since SNOOPY solves this set 
of equations. The shearing box equations in the frame co¬ 
moving with the background shear velocity v s * = -qClxe y are: 

tt + v s h -r~ + V ■(*'*' + T) = 2£lv y e x + (2 - q)Q\\e y + vV 2 v, (5) 
at dy 

— = V x (v x b) + qW 2 b, (6) 

at 

V ■ v = 0, (7) 

V ■ b = 0, (8) 

where v and b are the velocity and magnetic field respec¬ 
tively. Here T is a stress tensor given by 

T = {p + b 2 /2)I-bb, (9) 


where I is the identity matrix and p is thermal pressure. 

Upon volume averaging the Navier-Stokes equation (eq. 
5), we obtain two coupled equations for the volume averaged 
velocities (v x ) and (Vj,): 


d(v.r) 

dt 

d(v y ) 

dt 


= 2Cl(v y ), 

= &{v x )(q - 2). 


( 10 ) 

( 11 ) 


which yields the solution that both averaged velocities are 
proportional to exp(±z'At) for q < 2, or ~ exp(±/cf) for q > 2 
where k 2 = 2Q 2 (2 - q). The above analysis shows that the ‘x’ 
and ‘y’ mean velocities will grow exponentially, if perturbed, 
in the Rayleigh unstable regime q > 2. This growth is a phys¬ 
ical effect for finite perturbations. However if we set initial 
mean velocities to be zero the physical velocities should re¬ 
main such, but in simulations they can grow because of trun¬ 
cation errors. We verified this with the finite volume code 
ATHENA. The truncation errors seeded the mean velocities 
and they grew exponentially bringing the simulation to a 
halt in just a few shear times (l/fgfl)). 

We therefore chose to use the publicly available incom¬ 
pressible pseudospectral code SNOOPY, which has the im¬ 
portant property that the box averaged mean velocities do 
not grow throughout the duration of the simulation. This 
is because the nonlinear terms in the code are of the form 
(ik-v)v, and do not contribute when k = 0. Linear terms can 
only contribute to k = 0 mode evolution if the initial value 
for the fields at k = 0 is not set to zero, but we started all of 
our simulations without perturbations in this mode. 


3 NUMERICAL RESULTS 
3.1 Setup 

Using SNOOPY, we solve the incompressible hydrodynamic 
and MHD equations in the shearing box approximation. 
We solve the equations where the background shear has 
been subtracted out. SNOOPY utilizes the the 2/3 antialias¬ 
ing rule (Canuto et al. 2006). Shear periodic boundaries 
are remapped every f rem a P = L y l(qilL x ) (Umurhan & Regev 
2004). We define the Reynolds and magnetic Reynolds num¬ 
bers Re = LiqSl/v, Rm = LcqCl/q, respectively, where L z = Q = 
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Figure 1. Time history plot of kinetic energy (solid) and 
Reynolds stress (dotted) for hydl5 (q = 1.5, red), hyd21 (q = 2.1, 
green) and hyd42 (q = 4.2, blue). The y-axis is in log scale and 
the x-axis is time in units of I /Q. 

1 in code units. We fixed Re = Rm = 1600 for most of our 
runs. 

We use large scale noise as initial perturbations (with 
zero mean) and set the net initial vertical field Bg = 0.025 
in code units, which corresponds to an initial plasma beta 
13 = RQ. 2 /(Bq/ 2) = 3200. The magnetic field is calculated in 
Alfven speed units. For all of our runs, we use the domain 
size L x = Ly = L z = 1 with a resolution of 64 J . Table 1 provides 
a summary of our runs. 

3.2 Hydrodynamic shear flow stability 

As discussed in the previous section, the q < 2 regime is 
stable in hydrodynamics (see also Tillmark & Alfredsson 
(1992), Bech & Andersson (1997), Brethouwer (2005) for 
earlier work). We checked this by simulating the Keplerian 
q = 1.5 regime as well as two different values of shear in 
the Rayleigh unstable regime q = 2.1,4.2. We plot the time 
history of the kinetic energy and the Reynolds stresses in 
Fig. 1. As predicted by the standard modal analysis, the 
Keplerian flow is stable and its fluctuations exponentially 
decay to zero whereas the two Rayleigh unstable runs reach 
a saturated turbulent state in just a few shear times. 

3.3 MHD shear flow stability 

For MHD the regime 0 < q < 2 is unstable to the MRI. In 
Nauman & Blackman (2015), we focused on the dependence 
on q for q < 2 and found that the results were consistent 
with the linear calculations of Pessah et al. (2006) and the 
empirical results of Abramowicz et al. (1996). In contrast, 
the q > 2 case is stable to the MRI so a comparison of 
saturated states of the two regimes is instructive. 

One common feature visible from Figs. 1, 2 and 3 is that 
the case of largest shear (blue line, q = 4.2) has the largest 
growth rate in both magnetic and kinetic energies. The trend 
of increased growth rate with shear is also a property of 
the q < 2 (k > 0) MRI regime (Nauman & Blackman 2015). 
However, the important difference to note both in Fig. 2 and 
3 is that the growth rate of the kinetic energy (Reynolds 
stress) is greater than that of magnetic energy (Maxwell 
stress) in the q > 2 regime. 


Figure 2. Time history plot of kinetic (solid) and magnetic en¬ 
ergies (dotted) for mhdl5 (q = 1.5, red), mhd21 (q = 2.1, green) 
and mhd42 (q = 4.2, blue). 
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Figure 3. Same as Fig. 2 but for Reynolds (solid) and Maxwell 
stresses (dotted). 

To further explore the difference between kinetic and 
magnetic energy in the q > 2 regime, we increased Re and 
Rm to 6400 and 12800 (at Pr M = Rm/Re = 1) for q = 4.2 
and observed that the ratio of kinetic energy to magnetic 
energy in the saturated state decreased to nearly 2.7 for 
Re = Rm = 12800 compared to ~ 5.0 for the Re = Rm = 1600 
and 6400 cases. An extensive study of Re, Rm dependence is 
beyond the scope of the current paper. For Keplerian flow, 
the turbulent stresses also depend on dissipation coefficients 
(see for example, Riols et al. (2015)). 

As reviewed in section 2.2 above linear theory suggests 
that we can break the dispersion relation into two differ¬ 
ent types of modes (Shakura & Postnov 2015): Rayleigh and 
Velikhov-Chandrasekhar (VC). For q > 2, the VC mode is 
stable at all wave numbers. Our results show that for q < 2, 
the magnetic energy leads the kinetic energy while for q > 2 
the kinetic energy leads the magnetic energy. This result is 
reminiscent of isotropically forced box simulations of MHD 
turbulence in the following sense. In such simulations, the 
turbulent driver is imposed by hand as a forcing function. 
Normally the forcing is in the the Navier-Stokes equation, 
but it can also be imposed in the induction equation. When 
the forcing is imposed in the Navier-Stokes equation the sat¬ 
urated state reveals that the kinetic energy dominates the 
magnetic energy at the forcing scale and below. In contrast, 
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Run 

Shear 

C x v y ) 

-{Eby) 

Q'kin.v = <v*v y )/<v^> 

^magj — (bxby)/(b y) 

mhd 15 

1.5 

0.4837 ±0.3812 

3.3940 ± 3.8068 

0.4027 ±0.1751 

0.1877 ±0.0257 

mhd21 

2.1 

0.5896 ±0.5397 

0.4824 ±0.3851 

0.8701 ±0.4520 

0.2334 ± 0.0889 

mhd42 

4.2 

2.3140 ±3.3113 

0.0735 ±0.0871 

0.5987 ±0.1660 

0.1712 ± 0.1164 

hyd21 

2.1 

0.0519 ±0.0624 


0.7331 ±0.2993 


hyd42 

4.2 

0.9823 ± 0.9622 


0.5210 ± 0.1351 



Table 1. The first three runs are MRI runs whereas the last two are the purely hydrodynamic runs. We do not list the Keplerian 
hydrodynamic run here as it did not become turbulent. All the quantities are time averaged from 1000(l/£2) to 2000(1/12) (time averaging 
is defined by an overline) for all of the runs and volume averaged (represented by angled brackets) over the whole box. The stresses {v x v y ) 
and (b x b y ) are normalized by L?£2 2 , which equals unity according to our definitions. The fifth column represents the ratio of the Reynolds 
stress to the square of the azimuthal velocity Q'kin.y = (v x v y )/{v y ), while the last column shows this ratio corresponding to the magnetic 
field a ma g,y = ~{b x b y )j(b},). It appears that Q'kin.y is a sensitive function of the shear parameter while ar ma g,y is roughly constant. 


when the forcing is in the induction equation, the mag¬ 
netic energy dominates the kinetic energy at these scales 
(Park & Blackman 2012). 

These circumstances reflect the fact that the transfer 
of energy from the quantity that is driven (v or b) is not 
100% efficient to the response quantity (b and v, respec¬ 
tively). Interpreted in this way, the results from our simu¬ 
lations suggest that the for the q > 2 regime, the Rayleigh 
mode acts more like an an effective “driving” in the Navier 
Stokes equation, whereas for the q < 2 regime, the VC 
mode perhaps leads to a kind of “effective” forcing in the 
induction equation. This physical distinction may be use¬ 
ful in the path toward constructing analytic theoretical ap¬ 
proaches and is consistent with toy models in the MRI 
context that invoke forcing in the induction equation (e.g. 
Squire & Bhattacharjee (2015)). More work is needed to as¬ 
sess this rigorously. 

Finally, we note that boxes that are sufficiently large 
in the direction normal to the shear (L y ,L z » L x ) can lead 
to qualitatively different regime of ‘spatiotemporal chaos’ 
(Pomeau (1986),Philip & Manneville (2011)). For q < 2 
MHD shearing box simulations with L : » L x . Shi et al. 
(2016) showed that coherent structures appear in the mag¬ 
netic field while more recently Nauman & Pessah (2016) 
have shown that both velocity and magnetic fields develop 
coherent structures. The boxes used in the present study 
have L x = L y = L z = 1, so the extent to which a similar role 
of large boxes might also apply to the Rayleigh unstable 
regime should be investigated in future work. 


3.4 Correlation in space (x-y plane) 

Studying the physical effect of shear on the flow is aided by 
computing the autocorrelation function (ACF) of the veloc¬ 
ity and magnetic fields in the x-y plane. This autocorrelation 
provides a dimensionless of measure of the length or time 
scale over which the velocity (or magnetic field) maintains 
a value similar to itself and thus provides a measure of the 
locality of interactions in a turbulent flow. For random func¬ 
tions, the ACF decays exponentially. A plot of the spatial 
ACF in the x-y plane characterizes the spatial anisotropy 
of the velocity and magnetic field fluctuations. 

Following the convention used by Guan et al. (2009) 
and Simon et al. (2012), we define the spatial ACF of the 


magnetic field component ‘i’ ( i = x,y, or z) as: 


ACF(6(c5x)) = 


2; f b, ,(x + (5x, t)bj(x, t)d?x 
f b 2 (x, t)cPx 


( 12 ) 


where b 2 = b\ + b 2 + b\. Note that ACF(ft) is normalized to 
its maximum value at zero displacement (Sx = Sy = Sz = 0). 
Like Guan et al. (2009), we subtract off volume averaged 
mean quantities (b = b to ta.i - (b )). The overline represents the 
time averaging over ~ 1000(1/12) time units of the saturated 
state. We use the analogous definition for the autocorrelation 
of velocity fields ACF(v(<5x)). 

Fig. 4 shows the ACF(v(5x)) and ACF(/>((5x)) of the three 
shear values we study in this paper, q = 1.5,2.1,4.2 for both 
the hydrodynamic and the magnetohydrodynamic runs. In 
contrast to previous work on the MRI (e.g., Guan et al. 
(2009), Simon et al. (2012), Nauman & Blackman (2015)), 
the tilt angle observed in plots of ACF(b(Sx)) with respect 
to the y-axis is not constant with respect to variations in 
q. In addition, the hydrodynamic velocity ACF in fig. 4 for 
both q = 2.1,4.2 is more localized compared to the MHD 
counterparts at these same q. Comparing the MHD ACF 
plots, the q = 2.1 and 4.2 MHD runs show a very localized 
magnetic field compared to the q = 1.5 run. 

The tilt angles for the q < 2 cases previously studied 
were successfully modeled using an analysis of shear on fluc¬ 
tuations which assumed linear terms dominated nonlinear 
terms in the Navier-Stokes equation. Given that the q > 2 
cases studied here do not show the same simple monatonic 
dependencies, we are led to investigate how the ratio of non¬ 
linear to linear terms in the MHD equations vary a function 
of q. In the next section, we will show the non-linear terms 
in the Navier-Stokes equation do indeed dominate the linear 
terms for the q > 2 case when compared to the q <2 MRI 
unstable cases of previous work. This is a step toward iden¬ 
tifying the source of the more subtle dependence of tilt and 
localization on q in the q > 2 regime even if though exact 
dependence cannot yet be predicted analytically. 


3.5 Shear dependence of stress and energy: 

nonlinearities are more influential for q > 2 
than q <2 

Here we provide three lines of evidence consistent with non¬ 
linear terms being more influential than linear terms when 
it comes to understanding the behavior of stress and energy 
in saturation as a function of q for the q > 2 regime com- 
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Figure 4. Contour plots of the autocorrelation of velocity and magnetic fields for different runs. 


pared to the q < 2 regime. This is why it is more difficult 
to explain the q trends of tilt angle and localization for the 
q > 2 regime than the q <2 regime. 


3.5.1 Navier Stokes equation: Explicit comparison of 
nonlinear vs. linear terms for different q regimes 

To investigate the effect of shear on the turbulent properties 
of the flow, we study the time history of the energies and 
stresses at early times before the flow reaches nonlinear sat¬ 
uration. We focus on the the ‘x’ and ‘y’ velocity equations 
here: 

d,v x = 2£lv y + B 0 d : b x + vV 2 Vj, + b ■ Vb x - v • Vv* (13) 

d t Vy = (q - 2)Slv x + B 0 d-by + vV 2 v v + b ■ Wb y - v • Vy,. (14) 

The last two terms represent non-linear terms in both 
equations. For q = 2, eq. 14 has no source term in the lin¬ 
ear regime and is similar to the (non-rotating) plane Couette 
flow but with v x taking the role of shear velocity. In contrast, 


for q = 4 the source terms in eqs. 13 and 14 are both pro¬ 
portional to 2Q. The q = 4 case results in apparent isotropy 
in the two components for the linear regime. 

We plot the evolution of the different linear terms in the 
two equations and compare them with the rms value of the 
non-linear terms v • Vv and b ■ Wb for early times first 20fT 1 
times in figs. 5, 6, 7. 

For q = 1.5, the 2flv y term in eq. 13 is comparable to 
the non-linear terms for q = 1.5 (top left panel of fig. 5), 
suggesting that for q <2 the linear effects are very influential 
even as the saturated state is approached. This is assessed 
visually by noting that the red dashed curve overshoots the 
magnetic curve at most in the last few time steps of this 
plot. The linear term due to magnetic tension, B 0 d z b x i , is 


3 For an initially zero net flux case, such a term would be absent 
in the linear limit. We did carry out zero net flux simulations for 
Hz,ini = Bosink x x for all three shear values at Re = Rm = 1600 and 
found that only the q = 4.2 run shows growth and sustenance of 


MNRAS 000, 1-10 (2015) 


























6 F. Nauman and E. G. Blackman 


q=1.5 



q=1.5 



Time(£2 ’) 

Figure 5. The comparison of different linear and non-linear terms 
in eq. 13 (top panel) and 14 (bottom panel) for the first 50 
times, with q = 1.5. 

nearly an order of magnitude weaker than the other terms 
in this plot. 

For q > 2 the top panels of (figs. 6, 7), show that the 
corresponding linear terms are nearly an order of magni¬ 
tude weaker than nonlinear terms 13. Note here that the red 
dashed curve dominates over a longer range of time com¬ 
pared to the q = 1.5 case. Since the non-linear effects are 
dominating the linear velocity and magnetic field terms in 
this regime, the flow in this regime is expected to be more 
random with a smaller correlation length, consistent with 
fig. 4. Note also that for q > 2 (particularly in the q = 4.2 
plot) the non-linear magnetic terms b ■ Vbj (where i = x or y) 
are considerably weaker than the corresponding non-linear 
velocity term v • Vv,- (red dashed), suggesting that magnetic 
effects are subdominant in both the linear and non-linear 
regimes for the q > 2 regime (eq. 13). 

Analogously, comparing the linear vs nonlinear terms of 
14 for q - 1.5 vs q > 2 we find that in this case the nonlinear 
terms dominate the linear terms in both regimes, but that 
the red dashed curves of the bottom panels of figs. 6 and 

kinetic and magnetic energy while for the other two runs, both 
kinetic and magnetic energy decay. 


q=2.1 



q=2.1 



Figure 6. The comparison of different linear and non-linear terms 
in eq. 13 (top panel) and 14 (bottom panel) for the first 50 
times, with q = 2.1. 

7 are more dominant over a longer time range than in the 
bottom panel of fig. 5. 

3.5.2 Induction equation: Explicit comparison of 

nonlinear vs. linear terms for different q regimes 

The induction equation has the form: 

d,b x = B 0 d z v x + qV 2 b x + b ■ Wv x - v ■ \b x (15) 

d,b y = -qdb x + Bod,b y + qV 2 b y + b ■ Vv v - v ■ Vb y . (16) 

The first two terms in the b x equation (eq. 15) and the 
first three terms in the b y equation (eq. 16) are linear. The 
terms of the form v • Wb and b ■ Vv are nonlinear because 
the velocity fields depend on the magnetic fields through the 
Navier Stokes equation (eqs. 13 and 14). When the magnetic 
fields are weak b 2 « v 2 , then these terms could be considered 
approximately linear. However, for all of the shear values 
considered in this paper, the magnetic and kinetic energy 
are comparable right from the beginning of the simulations 
so it appears that the last two terms in both eqs. 15 and 16 
are nonlinear. 

The bottom panel in figures 8, 9 and 10 show that the 
generation of the azimuthal field b y due to the shearing of 
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q=4.2 



q=4.2 



Figure 7. The comparison of different linear and non-linear terms 
in eq. 13 (top panel) and 14 (bottom panel) for the first 50 f2 _1 
times, with q = 4.2. 


q=1.5 



Time(£2 ') 


q=1.5 



Figure 8. The comparison of different linear and non-linear terms 
in eq. 15 (top panel) and 16 (bottom panel) for the first 50 
times, with q = 1.5. 


the radial field b x is very significant in the first few rotation 
times flF 1 but is nearly an order of magnitude weaker than 
the v-Vby term in the saturated regime. The other nonlinear 
term b ■ Vv v is slightly larger in magnitude for q = 2.1 and 
q = 4.2 compared to the qQ.b x terms in the saturation regime 
but the two terms are nearly equal for q = 1.5. This suggests 
that stretching is more important for field growth in the 
q > 2 regime than the q <2 regime. 


3.5.3 Dependence of stresses and correlation time on q 

To evaluate how a kin!y (= {v x v y )!{vj)) and a m ^ y (= 
~{b x by)l{bl)) vary with shear, we use the autocorrelation 
function of time to obtain the correlation time. Following 
our earlier work (Nauman & Blackman 2015): 


ACF(v,(*» = 


" J v y (x, t + dt)v y (x, t)dt ' 

/ v 2 y(x,t)dt 


(17) 


where the angle brackets represent volume averaging over 
all space. Time integration is done over several orbits in 
the turbulent saturated state. Similarly we can calculate the 
cross correlation in time of b x with b y and v x and v y . For 


example, for the velocities we have: 


CCFTv^f*)) 


' ' 

J v x {x, t + 6t)v y (x, t)dt 

yj f vl(x, t)dt f v*(x, t)dt 


(18) 


The correlation times, computed from an exponen¬ 
tial fits to the plot of the ACF or CCF as in Fig. 11, 
tells us the characteristic time scale over which turbulent 
quantities such as the velocity are correlated to themselves 
or other quantities. From MRI simulations with q < 2, 
Nauman & Blackman (2015) found that the correlation time 
between x and y components of the field r was roughly 
inversely proportional to the shear. There the stress ACF 
was calculated instead of the CCF but we checked that the 
CCF exhibits a similar l/q behavior in the q < 2 regime 
Nauman & Blackman (2016). 

The importance of the correlation time is that when 
linear stretching in the induction equation can be used to 
estimate the amplification of azimuthal fluctuations from 
radial fluctuations, the azimuthal field is amplified by shear 
during a correlation time with dominant term 


O-rnag.y = ~ \q^\T (19) 


If r ~ 1/qCl, O'mag.y is roughly constant with shear. Indeed 
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q=2.1 



q-2-i 



Time(D 1 

Figure 9. The comparison of different linear and non-linear terms 
in eq. 15 (top panel) and 16 (bottom panel) for the first 50 12 
times, with q = 2.1. 


q=4.2 



q=4.2 



Figure 10. The comparison of different linear and non-linear 
terms in eq. 15 (top panel) and 16 (bottom panel) for the first 50 
IT 1 times, with q = 4.2. 


for q <2 case, this was confirmed by the simulations. More¬ 
over, the correlation times for the three quantities v y , b x and 
b x b y were very similar (see figure 13 of Nauman & Blackman 
(2015)). Using a similar argument for velocity as in eq. 19, 
we would get 

Ukin.y = (v x v y )/{v 2 y ) ~ (|(<? - 2)0|l)-'. (20) 

We now assess whether the above two equations, which 
are rooted in linear analysis, are equally effective at explain¬ 
ing the trends found in the q > 2 cases. We focus on the CCF 
(which is more relevant that the ACF) for stresses. We find 
that O'mag.y is nearly constant just like the q <2 MRI regime 
(see table 1) owing to the l/q dependence of the correlation 
time for CCF ( b x b y {8t )) (fig. 11). However, orkin.y decreases 
both for the HD and MHD runs unlike the q <2 cases 4 Eq. 
(20) would require that for Ukm.y to decrease, r has to go 

4 Fig. 7 of Nauman & Blackman (2015) shows that ( v x v y )l{v 2 ) in¬ 
creases with shear. We did not plot the CCF ( v x v y (6t )) in that 
paper but we checked that the correlation time for the Reynolds 
stress also varies as l/q, which in a linear picture would explain 
the increase in the ratio of Reynolds stress to the kinetic energy 
as eq. 20 suggests with the assumption (v 2 ) ~ (v 2 ). 



5t(f2 _1 ) 

Figure 11. The CCF (b x b y (6t)) as defined in eq. 18 but only for 
MHD runs. The x-axis is in units of 1/H. The colour scheme is as 
follows: mhdl5 (red), mhd21 (green), mhd42 (blue). 
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Figure 12. The ACF(v y (£?)) as defined in eq. 17. The colour 
scheme is as follows: mhdl5 (red), mhd21 (green), mhd42 (blue), 
hyd21 (magenta), hyd42 (black). 



Figure 13. The CCF (v x v y (6t)) as defined in eq. 18. Colour scheme 
same as fig. 12. 



q 

Figure 14. Correlation time calculated from an exponential fit 
to the MHD simulation plots in figs. 12, 13, 11. The y-axis is in 
units of 1/n. 


down faster than \q - 2|~‘. To check this, we plot the ACF of 
v y in fig. 12, which shows a slight increase with shear. Then 
r would be predicted to decrease by a factor of more than 22 
as q varies from q = 2.1 to q = 4.2 if Eq. (20) were the whole 
story. But the CCF of v^v,, in fig. 13 shows only a factor of 3 
(for MHD) to 4 (for HD) times decrease with shear (see fig. 
14). The likely explanation for this discrepancy is that Eq. 
(20) does not capture the effect of nonlinear terms. Indeed 
the comparison of linear and non-linear terms in figs. 6, 7 
shows that non-linear terms are generally more important 
than the linear terms for q > 2 regime. 

3.5.4 Tilt angle dependence on q 

The tilt angle in ACF [b(6x)) (fig. 4) has been di¬ 
rectly connected to the ratio of Maxwell stress to mag¬ 
netic energy {-b x b y )/(b 2 ) in previous work on MRI (e.g. 
Nauman & Blackman (2015)). Here we modify the definition 
to compare the stress to just the y-component of magnetic 
field squared, a mag0 , = ~(b x b y )/{b 2 ) = tan0 tilt . 

For our Rayleigh unstable simulations, the tilt angle 
observed from the ACF (b(5x)) and the definition based on 
n magJ , 5 disagree, in contrast the MRI q < 2 cases. For ex¬ 
ample, for q = 2.1 the a mag0 , = 0.2334 which is equivalent to 
#tiit ~ 13.14° whereas for q = 4.2, the n magJ , = 0.1712 which 
is equivalent to #tilt ~ 9 .71° (fig. 4). A visual inspection of 
fig. 4 shows that the q = 4.2 tilt angle is nearly 45°. 

From our discussion in sections 3.5.1 and 3.5.2, we are 
led again to the conclusions that this is further evidence for 
the more dominant role of nonlinear terms in the Rayleigh 
unstable regime compared to the MRI unstable q <2 regime. 
This demonstrates the inadequacy of linear arguments to 
explain the correlation between b x and b y . 

At present we do not have a non-linear model to explain 
the observed behavior in either the spatial correlation (fig. 4) 
and the temporal correlation (fig. 14) but the identification 
that the nonlinear terms are essential is a step toward such. 
The importance of these nonlinear terms present a challenge 
for analytic explanations. 


4 CONCLUSIONS 

We have compared the turbulent saturation properties of 
Rayleigh unstable MHD shear flows with those of the more 
commonly studied MRI unstable but Rayleigh stable regime. 
Our results are summarized below: 

(i) The Rayleigh unstable regime (q > 2) generates tur¬ 
bulent velocity flows with or without magnetic fields. In the 
presence of magnetic fields, the fluid turbulence drives dy¬ 
namo amplification of the total magnetic energy. 

(ii) In this q > 2 regime, we find that magnetic energy and 
Maxwell stresses saturate at lower values than the kinetic 
energy fluctuations and associated Reynolds stresses. In this 
regime therefore, the magnetic field is “slaved” to the flow 
turbulence. This contrasts the MRI unstable regime (q < 2) 
in which the magnetic fluctuations and magnetic stresses 
dominate the kinetic energy fluctuations and stresses. 


5 We thank the referee for pointing this out. 
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(iii) The quantity a ma g, y remains roughly constant in 
q > 2 regime, which is the same as for the q < 2 
regime (Nauman & Blackman (2015)). The tilt angle in 
ACF(b(6x)), on the other hand, with respect to y-axis is not 
constant as q changes. This contrasts the behavior in the 
MRI regime where the tilt angle is constant with changing 
<?■ 

(iv) We found that the magnetic structures of the flow 
become more localized as we increase the shear from q = 1.5 
to 4.2. 

Our work on MHD turbulence in the Rayleigh unstable 
regime has shown qualitative differences in the way quanti¬ 
ties scale with q compared to the more well studied MRI un¬ 
stable regimes. While the dependencies on q for MRI regime 
seems to be captured by analytic explanations that invoke 
linear analysis, the same linear estimates do not work for the 
q > 2 cases. We have traced the source of these differences to 
the stronger influence of non-linear effects in the Rayleigh 
unstable regime. A physical and analytic understanding of 
these differences requires non-linear modeling of MHD shear 
turbulence in the two regimes, which is good opportunity for 
work beyond the present scope. 
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